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ABSTRACT 

Granular  soil  is  represented  in  this  work  by  2-D  random 
arrays  of  elastic,  rough,  quartz  spheres  using  the  "distinct 
element"  method.  The  3-D  computer  code  TRUBAL,  originally 
developed  by  Peter  Cundall,  has  been  modified  at  RPI  by  the 
introduction  of  a  general  solution  to  the  Hertz-Mindlin  contact 
problem.  This  has  been  achieved  by  attaching  a  subroutine  to 
the  original  code,  which  describes  the  nonlinear 
force-displacement  relationship  at  the  intergranular  contacts, 
by  means  of  plasticity  theory  and  kinematic  hardening.  The 
above  modified  program  (CONBAL-2),  developed  as  part  of  another 
RPI  project,  can  perform  2-D  simulations  and  has  already  been 
used  to  study  the  dynamic  small  strain  behavior  as  well  as  the 
large  strain  behavior  of  sand.  The  2-D  small  strain 
simulations  in  this  report  include  the  determination  of 
compressional  wave  velocities  for  isotropically  and 
anisotropically  loaded  random  arrays  of  spheres.  The  results 
are  compared  with  excellent  agreement  to  existing  analytical 
and  numerical  solutions,  as  well  as  to  experiments  performed  on 
the  large  triaxial  cubical  specimen  at  the  University  of  Texas. 
It  is  found  that  the  fabric  of  the  material  plays  a  very 
important  role  in  the  above  phenomena.  In  all  stages  of  the 
simulations,  the  observed  macroscopic  phenomena  are  related  to 
the  microscopic  phenomena  at  the  interparticle  contact  level. 


INTRODUCTION 

Since  the  early  1960 's,  a  large  number  of  experimental 
studies  has  been  performed  on  the  small  strain  behavior  of  sand 
aiming  at  providing  empirical  correlations  for  practical  use 
(Hardin  and  Richart,  1963,  Hardin  and  Black,  1966,  Drnevich  and 
Richart,  1970,  Dobry  et  al,  1982).  The  equations  relating  the 
shear  modulus,  G__,,,  to  the  void  ratio,  e,  and  the  effective 
confining  pressure,  "a  ,  proposed  by  Hardin  and  Richart  (1963), 
Seed  and  Idriss  (1970)  and  Seed  et  al  (1984)  are  especially 
important  to  this  report.  These  empirical  relationships  are 
based  on  the  assumption  that  the  soil  can  be  treated  as  an 
elastic,  isotropic  solid,  in  all  three  correlations: 


m 


G=  A  f (e)  a? 
max  o 


where  A  is  a  constant,  f(e)  is  an  experimentally  defined 

function  of  the  void  ratio,  e,  ~5q  =  1/3 (a  1  +  T^)  is  the 

mean  effective  stress,  and  p  is  an  experimentally  defined 

exponent  which  has  been  found  to  be  approximately  0.5.  The 

above  empirical  relations  assume  that  the  value  of  G  (and 

max 

thus,  the  wave  velocity,  V  =  /G__/p,  where  p  is  the  mass 

S  IUoX 

density  of  the  soil),  is  the  same  for  isotropically  and 
anisotropically  loaded  sand,  provided  that  the  effective  mean 
stress,  a  ,  is  the  same.  Furthermore,  all  correlations  assume 


that  for  the  anisotropic  loading  case,  G _ 

max 

change  in  direction. 


and  V  do  not 
s 


a 


These  assumptions  for  Gmax  in  sands  have  been  challenged 

by  the  experimental  results  obtained  by  Schmertmann  (1978), 

Roesler  (1979),  Knox  et  al  (1982),  and  Yu  and  Richart  (1984). 

It  was  established  by  these  researchers  that  both  the 

compress ional ,  V  ,  and  shear,  Vg,  wave  velocities  depend  on  the 

stresses  in  the  direction  of  wave  propagation  and  particle 

motion  polarization,  and  not  on  the  effective  mean  stress. 

Stokoe  et  al  (1980)  developed  at  the  University  of  Texas 

at  Austin  a  large  scale,  7  X  7  X  7  ft  cubical  triaxial  facility 

for  the  specific  purpose  of  measuring  V  and  V  in  dry  sand. 

^  * 

In  this  facility,  a  triaxial  state  of  stress,  ai  *  °2  *  °3 
in  dry  sand)  can  be  achieved.  All  tests  performed  to  date  at 
the  University  of  Texas  have  used  a  locally  found,  medium  to 
fine,  washed  mortar  sand  classified  as  SP,  with  effective  grain 
size,  D^^O.28  mm  and  a  uniformity  coefficient,  Cu,  of  1.7. 
The  sand  is  placed  by  the  raining  technique  and  tested  dry. 
The  values  of  the  principal  stresses  have  ranged  between  15  psi 
and  40  psi,  with  the  stress  ratio,  K  =  *  1  to  2.67. 

Kopperman  et  al  (1982)  used  this  facility  to  study  the  effect 
of  the  stresses  on  the  propagation  velocity  of  compress ional 
waves,  V  .  They  caused  P-waves  to  propagate  along  one  of  the 


principal  stresses  (o  ),  while  varying  the  others;  they 

cl 

concluded  that; 

V  =  L  (fps)  (2) 

P  a 


where  L  =*  220  to  290  and  m  a  0.20  to  0.24  are  experimentally 

defined  constants.  The  sensitivity  of  V  on  a  and  its 

P  a 

n  . 

msesitivity  to  variations  m  the  other  two  principal  stresses, 
and  oc  perpendicular  to  wave  propagation  is  illustrated  in 
Fig.  1. 

In  a  similar  way,  Knox  et  al  (1982)  studied  the 
propagation  of  shear  waves  under  anisotropic  conditions  and 
concluded  that: 


Vs  = 


F  a 


ma  mb  me 


(fps) 


(3) 


where  F  is  a  constant,  ma  a  mb  a  0.09  to  0.12,  and  me  *  0  to 
0.01.  Thus,  in  these  two  cases  the  experimental  findings  of 
Schmertmann  (1978)  and  Roesler  (1979)  were  confirmed,  and  it 
was  established  that  empirical  correlations,  such  as  Eq.  1,  do 
not  adequately  describe  the  behavior  of  anisotropically  loaded 
sand,  and  they  need  to  be  revised  and  upgraded.  In  most  cases 
sand  is  anisotropically  loaded,  and  more  than  2  elastic 
constants  may  be  needed.  For  example,  if  ^  = 


o 2  *  o the  soil 


will  behave  as  a  transversely  isotropic  medium  and  5  constants 
(moduli)  will  be  needed;  in  the  more  general  case  of  a-±  *  °2  * 
a the  sand  may  behave  as  an  orthotropic  continuum  and  at 
least  9  constants  may  be  required  ( Sokolnikof f ,  1956). 

This  contradiction  between  the  early  and  the  more  recent 
laboratory  data  indicates  that  the  constitutive  modelling 
phenomenon  of  soils  at  small  strains  is  more  complex  than 
originally  anticipated.  While  it  is  true  that  a  number  of 
phenomenological  continuum  models  of  the  small  strain  (elastic) 
response  of  sand  have  been  proposed  (Rowe,  1971,  Coon  and  Evans 
1971,  Lade  and  Nelson,  1988),  they  apply  only  to  isotropic 
sands.  On  the  other  hand,  experimental  studies  by  Oda  (1974), 
Oda  and  Konishi  (1974)  ,  Oda  et  al  (1983,  1985)  have  indicated 
the  potential  of  a  particulate  mechanics  approach  to  explain 
the  observed  behavior.  This  is  further  strengthened  by  the 
findings  of  earlier  analytical  work  by  Duffy  and  Mindlin 
(1957),  Deresiewicz  (1958a),  Duffy  (1959),  and  others. 
Consequently,  it  is  generally  accepted  that  at  small  strains 
the  behavior  of  a  granular  medium  is  governed  by  the  elasticity 
of  the  particles  (i.e.  the  force-displacement  nonlinearity  at 
the  interparticle  contacts),  while  in  the  very  large  strain 
range  the  geometrical  arrangement  and  the  kinematics  of  the 
particles  become  dominant. 


Although  a  number  of  analytical,  experimental  and 
numerical  particulate  mechanics  studies  have  been  performed  in 
the  past  on  regular  and  random  arrays  of  spheres,  their  scope 
has  been  limited  by  the  available  computer  capabilities  and  the 
lack  of  a  realistic  force-deformation  law  at  the  intergranular 
contacts.  In  the  last  few  years  the  wider  availability  of 
supercomputers  for  basic  engineering  research  and  the 
development,  at  RPI ,  of  a  general  solution  to  the  Hertz-Mindlin 
problem  at  the  contact  between  two  identical  elastic,  rough 
spheres  (Seridi  and  Dobry,  1984,  Dobry  et  al,  1988)  have 
significantly  improved  the  situation.  This  solution,  coded  as 
program  CONTACT,  could  now  be  used  as  a  subroutine  in  a  finite 
element  program  which  describes  the  stress-strain  response  of 
granular  soil.  This  was  done  as  part  of  the  current  AFOSR 
sponsored  research,  and  Petrakis  and  Dobry  (1987,  1988a,  1988b) 
used  the  above  relationship  in  a  nonlinear  finite  element 
scheme  to  succesfully  predict  the  stress-strain  response  of  a 
regular/random  array  of  equal  spheres. 


RANDOM  ARRAYS  AND  TFE  DISTINCT  ELEMENT  METHOD 


Several  approaches  have  been  used  to  model  the  effect  of 
randomness  and  deviations  from  regularity  in  arrays  of  equal  or 
unequal  spheres.  Serrano  and  Rodrigues-Ortiz  (1973)  suggested 
a  method  for  generating  a  random  configuration  of  unequal  disks 
or  spheres  having  a  prescribed  grain  size  distribution. 
Cundall  and  Strack  (1979)  used  a  similar  approach  in  conjuction 
with  their  "distinct  element  method"  to  succesfully  simulate 
the  mechanical  response  of  arrays  and  disks  and  spheres  under  a 
variety  of  loading  conditions.  In  their  method,  an  explicit 
finite  difference  formulation  is  used  to  determine  the  static 
response  of  the  array  to  applied  strains  (program  TRUBAL)  or  to 
a  boundary  displacement  (program  BALL). 

TRUBAL  uses  a  periodic  space  in  a  form  of  a  cube  (3-D)  or 
square  (2-D),  such  as  that  of  Fig.  2,  to  minimize  the  effect  of 
the  boundaries  and  allow  the  use  of  a  relatively  small  number 
of  particles.  A  uniform  strain  field  is  applied  to  all 
particles  and  the  corresponding  contact  forces  are  calculated 
from  the  relative  displacements  between  neighboring  spheres. 
Initially,  the  program  used  an  arbitrary,  linear,  non  pressure 
dependent  law  to  describe  the  force-deformation  behavior  at 
each  interparticle  cointact;  in  a  later  version  of  the  program 
(Zhang  and  Cundall,  1986)  a  linear  pressure  dependent 


force-deformation  relationship  was  introduced.  As  illustrated 
by  Fig.  3,  these  contact  forces,  P^,  produce  in  each  sphere  a 
resultant  unbalanced  force,  EP^  ,  and  unbalanced  moment,  ZKL, 
which  induce  linear  and  angular  accelerations  in  the  particle. 
These  accelerations  are  used  in  turn  to  compute  a  new  particle 
position  at  the  end  of  a  time  increment,  At,  and  new  contact 
forces.  The  process  is  repeated  until  static  equilibrium  is 
achieved  (EKT  =*  a  0).  The  method  has  the  advantage  of 
decreasing  substantially  the  required  computer  memory,  as  no 
large  stiffness  matrices  need  to  be  calculated  and  inverted. 
However,  the  execution  time  is  large  due  to  the  great  number  of 
iterations  needed  to  assure  static  equilibrium. 

Cundall  and  Strack  (1979),  Strack  and  Cundall  (1984)  and 
Zhang  and  Cundall  (1986)  have  used  BALL  and  TRUBAL  for 
succesful  simulations  of  monotonic  direct  simple  shear  and 
compression  triaxial  tests. 

The  available  program  TRUBAL  (Strack  and  Cundall,  1984) 
was  modified  at  RPI  by  Ng  and  Dobry  (1988)  within  another,  NSF 
sponsored  project,  by  replacing  the  existing  arbitrary  linearly 
elast ic-perf ectly  plastic,  non  pressure  dependent, 
force-displacement  relation  at  the  intergranular  contacts,  by 
the  more  realistic  and  rigorous  Hertz-Mindlin  contact  law. 
This  was  accomplished  by  attaching  the  program  CONTACT 


described  previously,  as  a  subroutine  in  TRUBAL  (Table  1).  The 
new  modified  program  was  named  CONBAL-2  (CONtact  truBAL  in 
2-dimensions ) .  CONBAL-2  is  a  two-dimensional  version  of  TRUBAL 
which  incorporates  CONTACT  and  has  also  been  modified  to  run 
more  efficiently  on  the  computer.  Furthermore,  CONBAL-2  does 
not  allow  any  rotation  or  rolling  of  the  particles,  as  this 
causes  as  yet  unsolved  problems  in  the  numerical  simulations. 
The  2-D  character  and  lack  of  particle  rotation  in  CONBAL-2 
makes  the  array  somewhat  stronger  than  actual,  3-D  random 
arrays  and  soils;  Ng  and  Dobry  (1988)  have  approximately 
corrected  for  this  by  reducing  the  interparticle  angle  of 
friction.  Therefore,  CONBAL-2  is  very  similar  to  the  original 
program  TRUBAL,  except  for  the  two  differences  just  noted,  and 
for  the  very  important  introduction  of  the  effect  of  the  normal 
force  which  is  now  rigorously  modelled  by  the 
Hertz-Mindlin-Dobry  algorithm. 

The  accuracy  of  CONBAL-2  was  checked  by  comparing  the 
results  of  a  simulation  of  monotonic  pure  shear  loading  on  an 
isotropically  compressed  simple  cubic  array  of  quartz  spheres, 
to  the  rigorous  solution  obtained  by  Deresiewicz  (1958).  The 
"window"  of  the  nine  identical  spheres  used  in  the  CONBAL-2 
simulation  appears  in  Figure  4.  The  results  of  the  numerical 
simulation  are  plotted  in  Figure  4  (solid  line),  and  are  in 


excellent  agreement  with  the  rigorous  solution  (crosses). 

Additional  numerical  simulations  appear  in  Ng  and  Dobry  (1988); 

they  include  a  number  of  large  strain  numerical  experiments  on 

random  arrays  of  equal  or  unequal  quartz  spheres  under  cyclic 

and  static,  drained  and  undrained  conditions.  The  results  of 

these  simulations  are  in  excellent  agreement  with  comparable 

experimental  data  on  granular  soils  reported  in  the  literature. 

Moreover,  CONBAL-2  has  shown  to  have  considerable  predictive 

power.  For  example,  Ng  and  Dobry  (1988)  subjected  the  2-D 

array  of  57  unequal  spheres  shown  in  Fig.  5,  to  a  monotonical 

constant-volume  uniaxial  compression,  after  the  array  had  been 

_  2 

consolidated  under  an  isotropic  pressure,  aQ,  of  3.4  kgf/cm  . 
This  was  done  by  increasing  the  vertical  stress,  a^,  acting  on 
the  array,  while  varying  the  horizontal  compressive  stress,  a 2, 
as  much  as  necessary  to  keep  the  volume  constant.  To  model  the 
situation  during  the  corresponding  undrained  compression 
"triaxial"  test  in  which  the  "cell  pressure"  would  be  kept 
constant  at  oQ=3.4  kgf/cm  ,  a  "pore  water  pressure"  u=-(o2-3.4) 
was  defined. 

The  results  of  this  simulation  are  summarized  in  Fig.  6. 
The  figure  includes  the  calculated  axial  stress-strain  curve  in 
Fig.  6a,  the  variation  of  pore  pressure,  u,  vs  axial  strain  in 
Fig.  6b,  and  the  effective  stress  path  in  Fig.  6c.  The 
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behavior  of  the  granular  soil  simulated  is  clearly  dilative;  it 
is  accurately  modelled  by  CONBAL-2,  which  yields  results  very 
similar  to  those  measured  in  undrained  triaxial  tests  on 
uniform,  rounded,  dense  Ottawa  sand,  such  as  those  of  Vaid  and 
Chern  (1983)  (specimen  S-3),  also  included  in  Fig.  6.  At  first 
glance  they  may  not  look  comparable;  however,  the  scales  are 
different  and,  with  this  in  mind,  it  may  be  seen  that  the 
stress-strain,  pore  pressure-strain  curves  and  the  stress  path 
on  the  p-q  space,  are  remarkably  similar  between  specimen  S-3 
and  the  numerical  simulation. 


NUMERICAL  SIMULATIONS  OF  THE  SMALL  STRAIN  BEHAVIOR  OF  QUARTZ 
SAND 

After  program  CONBAL-2  was  checked  and  became  ready  to 
use,  a  copy  was  made  available  to  this  AFOSR  project  for 
research  on  the  stress-strain  behavior  and  development  of  a 
constitutive  relation  for  dry  granual  soil. 

As  a  first  step,  it  was  decided  to  start  simulating 
phenomena  in  the  small  strain  range,  and  at  a  later  stage  focus 
on  the  large  strain,  fully  nonlinear  inelastic  behavior. 
Furthermore,  it  was  also  decided  to  compare  the  results  of 
CONBAL-2  with  answers  obtained  using  the  analytical  and 
numerical  procedures  previously  developed  in  this  project  by 
Petrakis  and  Dobry  (1986,  1987),  as  well  as  with  experimental 
data  in  the  literature. 

Random  arrays  of  equal  spheres  were  generated  first,  since 
the  previous  analytical  and  numerical  methods  were  developed 
for  equal  spheres  and  much  available  experimental  data  are 
obtained  on  uniform  sand.  Moreover,  the  force-deformation 
relations  at  the  interparticle  contacts  developed  by  Hertz 
(1882),  Cattaneo  (1938),  Mindlin  (1949)  and  Mindlin  and 
Deresiewicz  (1953)  apply  only  to  equal  spheres  of  the  same 
material;  subroutine  CONTACT  is  based  on  these  formulations  and 
cannot,  in  principle,  model  the  contact  behavior  of  unequal  or 
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dissimilar  spheres.  Therefore,  program  CONBAL-2  yields 
rigorous  results  only  when  arrays  of  equal  spheres  are 
generated;  however,  the  same  program,  with  certain  assumptions, 
can  handle  in  first  approximation  moderately  unequal  spheres  of 
the  same  material,  provided  that  their  diameters  are  not  very 
different.  In  this  section,  results  on  packings  of  equal 
spheres  will  be  presented,  together  with  some  results  on  an 
array  of  moderately  unequal  spheres. 

i) Isotropically  Compressed  Random  Arrays  of  Spheres 

Two  different  2-D  random  arrays  of  477  identical  particles 

2 

having  the  properties  of  quartz  (G  »  295182  kgf/cm  ,  v  =0.15, 

f=0.5,  White,  1964),  were  generated  using  CONBAL-2.  In  the 

following  distinct  element  simulations,  f=0.35  was  used  instead 

of  f=0.5  to  correct  for  the  effects  of  2-D  and  absence  of 

particle  rotation  as  previously  discussed.  The  first  array 

was  a  very  loose  one,  with  an  average  coordination  number,  CN 

(number  of  contacts  per  sphere),  of  2.1,  while  the  second  was 

medium  dense  with  CN=3.0.  It  should  be  noted  that  when 

computing  the  average  number  of  contacts  per  sphere  in  a  random 

array,  all  spheres,  including  those  with  no  contacts,  are 

counted,  and  this  causes  the  average  coordination  number  to  be 

low.  These  two  arrays  were  then  subjected  to  three  different 

2 

values  of  isotropic  pressure  (aQ=  0.91,  3.34  and  6.98  kgf/cm  ), 


without  significant  change  in  their  average  number  of  contacts. 

The  configuration  of  the  second  assemblage,  with  CN*3.0, 

2 

and  subjected  to  °i1=ff22=  ao=  k9f/cm  is  shown  in  Fig.  7. 

The  rectangles  in  this  figure  represent  the  magnitudes  and 
directions  of  the  contact  forces;  there  are  four  different 
rectangle  widths,  each  one  of  them  corresponding  to  a  range  of 
forces  between  four  equal  fractions  of  the  maximum  computed 
contact  force.  For  example,  if  the  maximum  contact  force  is  F 
kgf,  the  narrowest  rectangle  stands  for  the  range  of  forces 
between  0  and  F/4  kgf,  the  next  wider  rectangle  for  the  range 
of  forces  between  F/4  and  F/2  kgf,  etc.  It  can  be  also 
observed  that  this  assemblage  of  equal  spheres  has  crystal ized; 
that  is,  areas  of  regular  packing  appear  within  the  random 
structure.  This  is  in  agreement  with  the  experimental  findings 
of  Smith  et  al  (1929),  Bernal  and  Mason  (1960),  Bernal  et  al 
(1964),  Davis  and  Deresiewicz  (1974),  Shahinpoor  and  Shahrpass 
(1982)  and  others,  reported  in  Petrakis  and  Dobry  (1986, 
1988a),  who  proposed  that  a  random  packing  of  equal  spheres  can 
be  represented  by  a  combination  of  regular  arrays.  Some  of 
these  regular  packings  in  Fig.  7  are  not  subjected  to  any 
contact  force,  with  the  forces  being  supported  through  arching 
by  surrounding  packings.  This  crystalizat ion  has  lumped  the 
array  into  N  constituents,  where  N  is  the  number  of  different 


packings  within  the  aggregate,  and  thus  has  increased  the 
characteristic  size  of  the  smallest  constituent.  As  a  result, 
the  array  is  not  completely  isotropic  under  isotropic  stress;  a 
much  larger  r  'mber  of  spheres  than  477  is  needed  in  order  to 
have  a  uniform  spatial  distribution  of  the  crystalized  regions 
and  achieve  a  statistically  isotropic  medium.  Moreover,  it  can 
be  seen  in  Pig.  7  that  the  applied  stresses  are  mainly 
transmitted  through  columns  of  particles,  in  which  all  contacts 
along  the  column  axis  experience  approximately  the  same  contact 
force.  These  contact  force  directions  cover  a  large  spectrum 
of  angles,  with  perhaps  some  preference  toward  the  two 
directions  parallel  to  the  boundaries.  These  conclusions 
obtained  from  a  visual  inspection  of  Fig.  7  are  also 
illustrated  in  Figure  8,  where  the  frequency  distribution  of 
contact  angles  is  shown  in  a  polar  plot  together  with  other 
micromechanical  statistics,  such  as  the  frequency  distributions 
of  mobilized  angle  (angle  between  contact  force  and  contact 
normal),  coordination  number  and  contact  force.  This 
477-particle  packing  can  be  approximated  as  isotropic  under  its 
current  isotropic  loading,  if  the  effect  of  the  orientation  of 
the  contacts  along  the  loading  directions  is  smoothed  in  Figure 
8. 
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The  two  477-particle  random  arrays  with  average 
coordination  numbers  CN=2.1  and  3.0,  were  subsequently  used  to 
compute  the  small  strain  shear  modulus,  G__.  under  isotropic 

luoX 

pressure.  For  this,  both  random  packings  were  subjected  to  a 

small  increment  of  macroscopic  shear  strain,  d7®lX10-^,  the 

corresponding  macroscopic  shear  stress,  dr,  was  calculated,  and 

the  shear  modulus  was  computed  as  G_  =dr/d7.  This  was  done 

for  both  arrays  at  the  three  confining  pressures  oq=  0.91,  3.34 
2 

and  6.98  kgf/cm  .  The  resulting  shear  modulus,  Gmax,  has  been 
plotted  versus  oq  and  the  coordination  number  in  Fig.  9,  which 
also  includes  the  small  strain  shear  modulus  for  the  simple 
cubic  array  in  two  dimensions  (CN*4).  For  a  given  oq  the  shear 
modulus  in  Figure  9  is  approximately  linearly  proportional  to 
CN,  similar  to  the  analytical  results  reported  by  Yanagisawa 
(1983)  and  Petrakis  and  Dobry  (1986,  1988a),  indicating  that 
the  stiffness  of  a  packing  of  spheres  is  controlled  by  both 
normal  stress  and  the  number  of  load  transmitting  contacts. 
Since  most  contacts  have  not  failed  in  the  random  array  of  Fig. 
7  (see  plot  of  mobilized  angle  in  Fig.  8)  the  whole  granular 
assembly  behaves,  through  its  contacts,  like  a  nonlinear  truss 
whose  stiffness  is  increased  when  new  members  are  added.  As 
expected,  the  deformation  characteristics  of  the  material  of 
the  particles  control  the  macroscopic  stress-strain  response  of 
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this  granular  medium  under  isotropic  pressure,  with  CN  being 
the  other  major  factor.  The  same  plot  may  also  be  interpreted 
as  a  plot  of  the  shear  wave  velocity,  V  ,  versus  the 
coordination  number,  CN,  if  it  is  assumed  that  the  wave  length 
is  significantly  larger  than  the  radius  of  the  spheres. 

Figure  10  contains  the  above  data  as  well  as  points 

derived  by  i)  the  Self  Consistent  Method  (Petrakis  and  Dobry, 

1986,  1988a);  ii)  analytical  solutions  for  regular  arrays;  and 

finally,  iii)  the  analytical  expressions  recently  reported  by 

Walton  (1937),  for  random  packings  of  equal  spheres.  Figure  10 

is  for  the  same  three  values  of  confining  pressure,  aQ  «  0.91, 

2 

3.34  and  6.98  kgf/cm  .  Walton,  by  considering  the  pressure 
dependent  force-deformation  behavior  of  both  normal  and 
tangential  compliances  at  the  interparticle  contacts,  derived 
the  following  expressions  for  the  two  elastic  moduli  (Lame 
constants),  X*  and  u*  -  Gmax'  of  a  random  packing  of  equal 
elastic  spheres  with  an  infinite  coefficient  of  friction, 


subjected  to  isotropic  compression,  oQ: 

*  C  3<^(CN)2<7o  1/3 

A  — - [ — - ] 

10(2B  +  C) 


*  (5B+C)  302(CN)2<7O  1/3 

A*  =  Gmax  =  “  ~~  [  JjTO  1 

10(  2B  +  C) 


(4) 

(5) 


**(l-n),  where  n  is  the  average  porosity  of  the  medium,  B  and  C 

are  constants  which  depend  on  the  material  properties  of  the 

spheres,  u  and  v  ,  and  CN  is  the  average  coordination  number. 

In  Fig.  10,  the  points  derived  by  the  Self-Consistent 

Method  and  Bqs.  4  and  5,  are  the  result  of  a  three  dimensional 

analysis,  unlike  the  points  in  Fig.  9  which  are  for  two 

dimensions.  Since  the  simple  cubic  array  appears  in  both  plots 

with  coordination  numbers  CN*4  (2-D)  and  CN*6  (3-D),  but  with 

* 

the  same  value  of  shear  modulus,  G  ,  it  was  decided  to 

luaX 

multiply  the  coordination  number  of  all  2-D  distinct  element 
simulations  by  6/4*1. 5,  to  approximately  "adjust"  for  the 
additional  third  dimension.  Furthermore,  since  the  3-D  Self 
Consistent  Method  was  applied  using  the  same  value  of  quartz 
modulus  for  the  spheres  as  used  in  the  other  methods,  but  with 
Poisson's  ratio  of  the  spheres,  v  *0,  in  order  to  eliminate  the 

S 

anisotropy  inherent  to  the  cubic  arrays  (Petrakis  and  Dobry, 

1986,  1988a),  Eqs  4  and  5  were  modified  accordingly. 

* 

Specifically,  for  »>  *0,  X  and  C  become  0  and  B*l/(2r«  ).  The 
final  result  is  the  very  consistent  plot  in  Fig.  10  using 
points  from  four  different  methods;  they  essentially  confirm 
the  hypothesis  presented  by  Yanagisawa  (1983)  and  Petrakis  and 
Dobry  (1986)  that  the  small  strain  shear  modulus,  G  ,  of  a 
random  array  of  equal  spheres  is  essentially  a  linear  function 


of  the  average  number  of  contacts  per  particle.  In  Eq.  5  the 

relation  between  G  and  CN  is  not  linear  but  the  product  of 

max 

the  multiplication  of  the  two  factors:  porosity,  n,  and  CN, 

* 

results  in  an  essentially  linear  relationship  between  u  =G  „ 

max 

and  CN,  once  the  influence  of  n  on  CN  is  considered. 

ii)Anisotropicallv  Compressed  Random  Arrays  of  Spheres 

The  same  array  already  discussed,  with  477  equal,  elastic 

spheres  and  coordination  number,  CN*3,  consolidated 

2 

isotropically  at  =  o22  =  0.91  kgf/cm  (Fig.  7)  was 

2 

further  loaded  under  biaxial  compression  to  a22  *  2.33  kgf/cm  , 
while  keeping  a ^  constant.  The  stress  path  appears  in  Fig,  11 
as  stress  path  (b).  This  was  done  to:  i)  investigate  the 
influence  of  the  magnitude  and  direction  of  the  principal 
stresses  on  the  velocity  of  longitudinal  waves,  V^,  propagating 
through  this  medium,  ii)  interpret  the  experimental  findings  of 
Kopperman  et  al  (1982),  and  iii)  verify  the  findings  of  the 
nonlinear  finite  element  model  presented  by  Petrakis  and  Dobry 
(1987,  1988b). 

To  compute  the  wave  velocities,  V  ,  the  constrained  moduli 
2 

D^-Vp  p,  were  calculated  first  as  follows:  once  the  desired 
stress  ratio,  ^■sai/'°2sa22/'all’  was  reac^e<^  at  specific  points 
along  stress  path  (b)  in  Fig.  11,  very  small  strain  increments, 
d«n,  d«22'  e<Iual  to  1.062X10-6  were  applied  to  the  array  with 
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appropriate  signs.  The  corresponding  stress  increments  do.^, 

d®22  wer®  computed,  and  finally,  the  small  strain  constrained 
(K) 

moduli  Dii  ”^ai i^^i i '  were  calculated  both  in  the  direction  of 

the  major  (°^*022^  an<^  minor  ^°2=an^  principal  stresses.  The 

results  of  these  simulations  are  shown  in  Fig.  12,  where  the 

small  strain  constrained  moduli  of  the  random  packing  at  every 

stress  ratio  K,  ,  normalized  by  the  constrained  modulus 

under  the  initial  isotropic  pressure  d| V ,  are  plotted  against 

the  stress  ratio,  K*a22  ^  all~°l/f°2 '  F*9ure  12  also  includes 

data  points  from  a  number  of  experiments  on  sand  performed  at 

the  University  of  Texas  by  Kopperman  et  al  (1982),  and  results 

of  the  nonlinear  finite  element  model  proposed  by  Petrakis  and 

Dobry  (1987),  all  performed  under  analogous  conditions  to  those 

of  the  CONBAL-2  nonlinear  distinct  element  simulation. 

The  agreement  between  all  normalized  results  in  Figure  12, 

is  excellent,  and  a  straight  line  of  equation  D22^D22^  = 

38  38 

^°22^all^*  *  K*  could  be  fitted  to  all  numerical  and 

experimental  data  points  included  in  the  figure.  It  should  be 
noted  that,  if  the  change  in  sand  density  as  K  increases  is 
neglected,  D22^D22^  =  [V^Vv^1^]2.  From  their  experiments  in 
sand,  Kopperman  et  al  (1982)  found  that  Eq.  2  was  applicable 
with  220<L*290  and  m*0.2.  This  gives  D^/D^^K0'4,  very 
similar  to  the  line  in  Fig.  12.  Furthermore,  it  was  found  that 
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the  actual  values  of  calculated  by  CONBAL-2  for  the  array  of 
Fig.  7,  are  within  10%  of  those  measured  by  Kopperman  et  al 
(1982).  Figure  12  shows  that  the  nain  conclusion  obtained  from 
the  University  of  Texas  laboratory  results,  that  the  P-wave 
velocity  (V  =  ✓D/p)  propagating  along  a  principal  stress 
direction  is  only  a  function  of  that  principal  stress,  is  fully 
predicted  by  the  nonlinear  distinct  element  method.  Again, 
this  is  an  effect  of  the  particulate  nature  of  the  soil,  which 
can  not  be  interpreted  and  reproduced  analytically  unless  this 
particulate  nature  is  taken  into  account. 

Figures  13  and  14  depict  the  positions  of  the  spheres,  the 

relative  magnitudes  and  directions  of  the  contact  forces,  and 

the  associated  micromechanical  statistical  information,  at  the 

last  stage  of  anisotropic  loading,  at  K=ff22//allIB2’'^//^*^1=^*^' 

This  is  the  last  stage  of  the  biaxial  loading  simulation,  whose 

particle  configuration  and  micromechanical  statistics  under  the 

2 

initial  isotropic  pressure  aQ=  an=a22=®’^  kgf/cm  (K=l)  was 
shown  in  Figs.  7  and  8. 

Comparing  the  array  configurations  of  the  above  two 
loading  stages  in  Figs.  7  and  13,  it  can  be  seen  that  the  array 
remains  crystal ized  throughout  the  loading  process,  and  that  a 
number  of  the  regular  packings  within  the  array  remain  stress 
free. 
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Similarly  to  Fig.  7,  the  four  different  widths  of  the 
rectangles  in  Fig.  13  represent  the  ranges  of  intensity  of 
contact  forces  relative  to  the  maximum  contact  force  in  each 
case.  When  the  random  array  is  subjected  to  isotropic 
pressure,  the  distribution  of  the  contact  forces  was  more  or 
less  uniform,  with  similar  contact  forces  in  all  directions 
(Fig.  8).  On  the  other  hand,  and  as  shown  in  Fig.  13  and  14, 
after  the  array  was  loaded  vertically  while  the  horizontal 
stress,  o^,  was  kept  constant,  the  number  of  contacts  in  the 
vertical  direction  increased  and  most  of  the  columns 
transmitting  loads  are  now  essentially  vertical  or  horizontal. 
These  columns  are  similar  to  the  "stiff  chains"  reported  by 
Cundall  and  Strack  (1983).  The  number  of  contacts  in  the 
horizontal  direction  has  somewhat  decreased,  while  in  all  other 
directions  has  decreased  substantially.  The  magnitudes  of  the 
contact  forces  in  the  vertical  contacts  have  also  increased. 

Figure  14  contains  the  micromechanical  statistics 
calculated  for  the  anisotropic  stress  condition  of  Fig.  13, 
including  frequency  distributions  of  contact  angle,  mobilized 
angle,  coordination  number  and  magnitude  of  contact  force.  The 
polar  plot  of  frequency  distribution  of  contact  angle, 
quantifies  the  increase  in  the  number  of  contacts  in  the 
vertical  direction,  previously  discussed,  as  well  as  the 
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decrease  of  number  of  contacts  in  the  other  directions.  Figure 
14  also  shows  that  at  this  point  of  the  biaxial  compression 
loading  most  contacts  have  failed,  the  particles  have 
rearranged  themselves  to  more  stable  conditions,  and  the 
distribution  of  coordination  numbers  has  changed  little.  The 
polar  plot  of  the  frequency  distribution  of  contact  angle  in 
Fig.  14  resembles  a  cross,  with  most  of  the  contacts  aligned 
parallel  to  the  applied  principal  stresses.  This  is  consistent 
with  the  vertical/horizontal  columns  transmitting  loads  in  Fig. 
13,  and  is  different  from  the  structure  of  the  array  under 
isotropic  stress  presented  in  Figs.  7  and  8.  There,  the 
frequency  distributions  of  contact  angle  and  contact  forces 
were  approximately  uniform  and  few  contacts  had  failed. 

Based  on  the  results  of  the  numerical  simulations  and  the 
above  discussion,  it  is  concluded  that  changes  in  the  structure 
of  the  random  packings  and  granular  soils  are  responsible  for 
the  observed  macroscopic  behavior,  and  specifically  for  the 
dependence  of  the  compressional  wave  velocity  only  on  the 
principal  stress  in  the  direction  of  wave  propagation.  Under 
isotropic  loading,  the  orientations  of  the  contacts  have  a  more 
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(3-D)  compression,  a  significant  number  of  contacts  is  gained 
in  the  direction  of  the  increased  major  principal  stress,  o^, 
while  contacts  are  lost  in  all  other  directions.  In  the  other 


two  principal  directions  in  which  the  stresses  a 2  (and  o^)  are 
kept  constant,  the  number  of  contacts  decreases,  but  not 
considerably,  while  many  contacts  are  lost  in  all  non  principal 
directions.  These  contacts  along  principal  directions  form 
"branches"  or  "stiff  chains"  parallel  to  the  principal  stresses 
which  transmit  most  of  the  applied  principal  stresses  from 
boundary  to  boundary.  The  result  is  a  structure  reminiscent  of 
a  simple  cubic  array  loaded  by  stresses  a^,  a 2  and  o 3  parallel 
to  the  main  axes  of  the  array,  and  with  the  effect  of  each 
principal  stress  being  uncoupled  from  the  other  two.  This 
uncoupling,  which  in  the  simple  cubic  array  happens  naturally 
due  to  th«»  geometry  of  the  array,  develops  in  the  random  array 
as  a  cc  equence  of  these  stiff  columns  or  "branches"  of 
particles  which  carry  the  applied  load.  Therefore,  as  in  the 
case  of  the  simple  cubic  array  (Petrakis  and  Dobry  1986,  1987), 
the  longitudinal  modulus,  D,  in  any  principal  direction  is  a 
function  only  of  the  stress  in  this  direction  and  is  unaffected 
by  variations  in  the  stresses  in  the  other  directions. 
Although  the  CONBAL-2  runs  discussed  here  were  two  dimensional, 
the  authors  postulate  that  the  same  phenomenon  occurs  in 


three-dimensions  with  the  formation  of  stiff  chains  in  three 
directions  instead  of  two.  It  is  well  known  (Deresiewicz  1958) 
that  the  shear  moduli  of  the  simple  cubic  array  depend  only  on 
two  principal  stresses  and  are  independent  of  the  third 
principal  stress.  Since  the  moduli  of  the  random  array  are 
expected  to  be  of  similar  nature  due  to  the  3-D  "stiff  columns" 
just  postulated,  it  is  believed  that  this  is  the  reason  why  the 
shear  moduli  (or  S-wave  velocities)  have  been  observed  to 
depend  only  on  the  principal  stresses  in  the  direction  of  wave 
propagation  and  particle  motion.  However,  this  hypothesis 
should  be  verified  once  the  3-D  version  of  CONBAL-2  is 
developed. 

To  further  demonstrate  the  fact  that  the  P-wave 
propagation  velocity  (o.r  the  constrained  modulus)  is  a  function 
of  the  principal  stress  in  that  direction  only,  the  same  random 
array  of  477  equal  spheres  was  subjected  to  a  biaxial 
compression-extension  loading  simulation  under  constant  mean 
stress.  oq»0,91  kgf/cm  .  The  corresponding  stress  path  (a)  is 
shown  in  Pig.  11.  The  vertical  stress,  °22’  was  ^ncrease<*  UP 
to  1.374,  while  the  horizontal  stress  was  decreased  down  to 
0.447  kgf/cm2  (Fig.  15).  Thus,  ctq=  1/2(1. 374+0. 447)*0. 91  as 
before.  Again,  the  normalized  constrained  moduli  d|^/d|^, 
were  computed  at  specific  values  of  the  stress  ratio,  K=o.^/ao, 
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DH^“D22^'  and  K=CTll^°o  for  °|i)“Dn)-  If  the  moduli  were  a 
function  of  the  mean  stress  (i.e  if  Eq.  1  were  true),  the 

P-wave  propagation  velocity  (and  the  modulus)  would  be 

unaffected  by  the  changes  in  the  magnitude  of  the  principal 

stresses.  would  be  always  the  same  in  both  directions, 

( K ) 

and  the  plot  of  D22  'D22  versus  K  *n  Fig*  16  would  be  a 
horizontal  straight  line.  The  normalized  constrained  moduli 
calculated  with  CONBAL-2  are  plotted  vs.  K  in  Fig.  16.  As  <722 
increases,  Dl-  increases  in  the  figure.  Simultaneously,  as 

decreases,  decreases,  and  the  plot  of  /D^  versus 

K  follows  a  consistent  curve  valid  for  both  directions. 
Therefore,  Fig.  16  confirms  again  that  the  constrained  moduli, 
are  affected  by  the  principal  stress  in  the  direction  of 
loading  only  and  are  completely  independent  of  the  mean  stress. 
Therefore,  based  on  the  above  results  and  those  of  the  biaxial 
loading  simulation  (Fig,  12),  it  is  predicted  that  a  granular 
soil  does  not  behave  as  a  linearly  elastic  isotropic  solid  at 
small  strains,  but  it  changes  its  order  of  elastic  symmetry 
depending  on  the  applied  principal  stresses.  Consequently, 
expressions  of  the  type  of  Eq.  1  can  be  used  only  for 
isotropically  consolidated  soils?  in  all  other  cases,  relations 


of  the  type  of  Eqs.  2  and  3  should  be  used.  The  contact  force 
distribution,  formation  of  stiff  columns  and  micromechanical 
statistics  of  this  array  at  the  end  of  the  loading  path  (Figs. 
15  and  17)  are  very  similar  to  those  obtained  in  biaxial 
compression  under  different  anisotropic  stress  conditions 
(Figure  13  and  14).  The  crystalization  appears  once  more; 
since  most  initial  contacts  have  failed  at  this  point,  new  have 
formed  along  the  direction  of  the  major  applied  principal 
stress.  At  the  same  time,  there  has  been  a  substantial 
decrease  in  the  number  of  contacts  in  the  horizontal  direction, 
and  stiff  columns  of  particles  have  appeared  both  vertically 
and  horizontally. 

Finally,  in  order  to  be  able  to  generalize  the  above 
findings  for  the  case  of  uniform  random  arrays  of  unequal 
quartz  spheres,  a  random  array  of  531  particles  of  two 
different  diameters  was  generated.  This  aggregate  was  composed 
of  168  spheres  of  radius  R1  and  of  363  spheres  of  I^.  The 
ratio  R^/Rj  was  set  to  1.5,  and  thus  not  very  different  from 
unity,  so  that  the  Mindlin-Deresiewicz  theory  could  be  applied 
by  using  the  concept  of  "equivalent  radius".  This  equivalent 
radius,  Rg^R-^/fR-^+Rj) >  where  R1  and  Rj  are  the  radii  of  the 
two  particles  in  contact,  has  been  used  in  the  past  in  the  case 
of  unequal  cylinders  (Poritsky,  1950),  and  represents  an 
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approximation  to  an  unknown  exact  solution.  It  is  derived  from 
the  Hertz  (1882)  theory  for  smooth  spheres  (f=0)  and  can  be 
applied  to  the  case  of  rough  spheres  as  a  first  approximation 
(Ng  and  Dobry,  1988) 

The  above  array  was  first  subjected  to  isotropic  pressure, 
2 

1.31  kgf/cm  ;  the  corresponding  geometric  and  statistical 
information  is  shown  in  Figs.  18  and  19.  As  the  spheres  are 
not  identical,  now  there  is  no  crystalization  and  the  531 
sphere  aggregate  is  almost  isotropic,  as  can  be  seen  from  the 
uniform  distribution  of  contact  angles  in  the  polar  plot  of 
Fig.  19  (compare  with  Fig.  8).  Therefore,  this  531-particle 
array  exhibits  a  higher  level  of  symmetry  than  the  477-sphere 
packing,  and  is  taken  to  be  isotropic.  Even  in  this  case,  in 
which  there  is  no  crystalization,  there  is  still  a  number  of 
particles  which  are  not  in  contact  and  do  not  carry  any  load. 
Moreover,  when  comparing  this  packing  in  Fig.  18  to  the  477 
equal  sphere  assemblage  in  Fig.  7,  no  clear  system  of  parallel 
columns  of  particles  can  be  detected  transfering  approximately 
the  same  forces  across  the  whole  array,  but  instead,  the 
magnitude  of  the  contact  force  varies  with  position. 

The  array  was  finally  loaded  under  biaxial  compression  to 
2 

022*3.32  kgf/cm  ,  while  keeping  o^  constant  as  shown  by  stress 
path  (c)  in  Fig.  11.  The  array  configuration  and  the  contact 
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forces  at  the  end  of  the  loading  path  are  shown  in  Figure  20, 


with  the  corresponding  micromechanical  statistics  presented  in 
Fig.  21.  As  we  can  see  in  both  figures,  the  contacts  have  been 
more  or  less  eliminated  in  all  directions  except  the  major  and 


minor  stress 

directions 

°l=ff22 

and 

in¬ 

The 

number  of 

contacts 

has 

doubled  in 

the 

ff22 

direction 

and 

has  also 

increased 

more 

moderately 

in  the 

°11 

direction. 

Most  contacts 

have  failed,  and  again,  there  are  still  some  spheres  which  are 
not  subjected  to  any  contact  force.  Finally,  the  contact 
forces  in  this  case  vary  again  with  the  position  of  each 
particle.  Vertical  and  horizontal  load  transmitting  columns  or 
"stiff  chains"  are  present  in  Fig.  20,  but  they  are  less 
continuous  and  not  so  obvious  as  they  were  in  the  crystalized 
equal  sphere  array  of  Fig.  13.  A  close  comparison  between  this 
array  of  unequal  spheres  in  Fig.  20  to  the  equal  sphere  array 
of  Fig.  13,  shows  that  while  their  structures  (fabrics)  are 
fundamentally  dissimilar  (crystalization  vs.  randomness)  and 
the  distribution  of  forces  within  them  very  different,  they 
both  form  the  majority  of  their  contacts  along  the  directions 
of  the  applied  principal  stresses,  when  compressed  biaxially. 
It  is  also  interesting  that  the  frequency  distributions  of  the 
contact  angle  in  all  simulations  at  the  end  of  their 
corresponding  biaxial  loading  simulations  (Figs.  14,  17,  and 


21)  have  the  distinctive  "cross"  pattern  with  the  contacts 
aligned  along  the  directions  of  the  applied  principal  stresses. 

The  normalized  constrained  moduli  were  also  computed  for 
the  array  of  Fig.  20  at  various  values  of  the  stress  ratio, 
anc*  they  are  Plotte£3  versus  K  in  Fig.  22.  Once  more 
the  moduli  (and  consequently,  the  P-wave  velocity)  are  affected 
in  first  approximation  by  the  principal  stress  in  that 
direction.  In  addition  to  the  results  of  the  simulation  on  the 
unequal  sphere  array  of  Fig.  20,  Fig.  22  also  includes  all 
points  of  Fig.  11,  that  is:  the  experimental  data  of  Professor 
Stokoe  and  his  co-workers,  the  nonlinear  finite  element  results 
of  Petrakis  and  Dobry  (1987,  1988),  and  finally  the  results  of 
the  nonlinear  distinct  element  method  on  the  array  of  the  477 
equal  eqraf^spheres.  All  points  are  in  excellent  agreement, 
thus  verifying  the  accuracy  and  versatility  of  the  approach, 
and  the  validity  of  the  conclusions  reached  for  the  constrained 
modulus,  D^,  and  P-wave  velocity,  V  ,  of  random  arrays  of 
spheres  and  of  granular  soils. 


CONCLUSION 


A  2-D  nonlinear  distinct  element  simulation  has  been 
presented  which  interprets  the  behavior  of  granular  soil  at 
very  small  strains.  This  simulation  is  based  on  an  incremental 
solution  to  the  nonlinear  problem  of  two  spheres  in  contact 
(program  CONTACT),  incorporated  into  the  available  distinct 
element  program  TRUBAL.  It  has  been  found  that  this  approach 
not  only  interprets  succesfully  the  small  strain  behavior  of 
granular  soil,  but  it  also  provides  a  wealth  of  information  on 
the  fabric  changes  during  loading  which  were  until  now  very 
difficult  to  obtain.  The  results  of  these  simulations  on  both 
isotropically  and  anisotropically  compressed  random  arrays  of 
equal  and  unequal  spheres,  are  in  excellent  agreement  with  a 
number  of  previous  analytical  and  numerical  procedures,  as  well 
as  with  experimental  data  on  sand  in  the  literature.  The 
origin  of  the  above  phenomena  is  the  changes  in  the  structure 
of  the  granular  medium  which  occur  when  loaded,  and 
specifically,  the  formation  of  columns  of  particles  supporting 
approximately  the  same  contact  force  along  the  directions  of 
the  applied  principal  stresses.  Therefore,  the  distribution 
and  magnitude  of  the  contact  forces  are  of  great  importance  to 
the  macroscopic  response  of  the  medium.  The  above  phenomena, 
which  are  micromechanical  in  nature,  can  only  be  modelled 
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TABLE  1 

(after  Ng  and  Dobry,  1988) 

RPI  COMPUTER  PROGRAM 


r  max. 
1500  I-  +1  SD 


Stress  Along  Axis  of  Propagation,  o^,  psi 


a)  North-South  Axis 


*  1800 


•  V500 


Stress  Perpendicular  to  Axis  of  Propagation,  oc>  psi 


b)  Vertical  Axis 


o  » 

>  .  1500  -  4— 

«  a 

5*  -£tfos- 

l  1300£ - 


Stress  Perpendicular  to  Axis  of  Propagation,  e»c,  psi 


c)  East-West  Axis 


Figure  1.  Effect  of  Structural  Anisotropy  on  Compression  Wave 
Velocity  for  Triaxial  Confinement  (after  Kopperman  et 
al,  1982). 
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DISTINCT  ELEMENT  METHOD 

(Programs  BALL,  TRUBAL  developed  by  Peter  Cundall) 


Iimz-1 


Pi 


Pi 


rZMi, 


2  Pi  produces  acceleration  of  particle 
ZMi  produces  angular  acceleration  (spin)  of 
particle 

Acceleration  and  Spin  Computed  as  if  Neighboring 
Particles  did  not  exist.  Time  step  is  assumed  to  be 
small  such  that  accelerations  and  velocities  are 
constant  in  this  small  time. 


•  Particle  occupies  new  position  due  to  acceleration 

•  New  Pi's,  Mi's  are  computed  for  new  positions 

•  rPi,  I  Mi  and  process  is  repeated. 


Figure  3.  Outline  of  the  Main  Features  of  the  Distinct 
Element  Method 
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Pigure  5.  A  Random,  Periodic,  2-D  Medium  of  57  Elastic,  Rough, 
Quartz,  Spheres,  of  two  Radii  (R./R^l.S),  Subjected  to 
Biaxial  Compression.  Note  tire  ^"window"  with  the 
representative  random  pattern  (after  Ng  and  Dobry, 
1988) 


s  •*  \  %  _ 


Pressure,  u,  Kgf/cm*  Axial  Stress 


Tirv0-91  k9  f/o» 


Figure  7.  2-D  Random  Array  of  477  Equal,  Elastic,  Rough,  Quartz 
Spheres  Subjected  to  Isotropic  Compression,  a  *0.91 
kgf/cm  .  Note  that  this  figure  represents  the  "window" 
with  the  representative  random  pattern. 
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Figure  8.  Statistical  Information,  Regarding  the  Isotropic 
Compression  at  o  -0.91  kgf/cnT  of  the  477-sphere  medium 
of  Figure  7.  ° 


Figure  9.  Shear  Modulus,  G  ,  Versus  Coordination  Number  (CN 
■  Number  of  Contacts’ prer  Sphere)  for  two  Random  Arrays 
of  477  Equal,  Quartz  Spheres,  and  the  Simple  Cubic 
Array  (CN«4) . 


Regular  Cubic  Arrays 
Nonlinear  Distinct  Element  Method 


Self  Consistent  Method  (Petrakis  &  Dobry,  1986) 
Walton  (1987) 


°o  =  6.98  kgf/cm 


on  =  3.34  kaf/cm 
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Figure  11.  Stress  Paths  Used  in  the  Simulations:  (a)Biaxial 
Compression-Extension  vith  Constant  Mean  Stress, 
o  *0.91  kgf/cnr;  (b)Biaxial  Compression  from  a  -0.91 
k§f/<-*-  — J  - - -  ' -  " 


(c)Biaxial  Compression,  from  o  *1.31 


Stress  Ratio,  K  *  c^^0!!  *  tV<J2 


1.0  1.5  2.0  2.5  3.0 


Stress  Ratio,  K  «  a22/,°n  “  ai/°2 


Figure  ’  Normalized  Constrained  Moduli,  D.  .  /D..  ,  Versus 
i.„.ess  Ratio,  K»o./a_  ■  a22^°ll'  Comparison  Between 

Distinct  Elevment,  Finite^  Element  and  experimental 
Results:  (a)  In  the  Direction  of  o2,,  and  (b)in  the 
Direction  of  o^  (o^  is  kept  constanrr. 
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Figure  14.  Statistical  Information  Regarding  the  Biaxial 
Loading  in  Figure  13. 


Q  in  the  direction  of 
*  in  the  direction  of  a ^ 


Figure  16.  Normalized  Constrained  Moduli,  D.,  /D..  ,  Versus 
the  Stress  Ratio,  K»c../o  ,  *or  the  ^®se  of  Biaxial 
Compression-Extension  .^oaaing  with  Constant  Mean 
Stress,  a  -0.91  kgf/cni 
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Figure  17.  Statistical  Information  for  the  of  Biaxial 
Compression-Extension  Loading  in  Figure  15. 


Figure  18.  2-D  Random  Array  of  531  Elastic,  Quartz  Spheres  of 
Two  Radii  (R./R--1.5)  ,  Subjected  to  Isotropic 
Compression  oQ*  IT 32*kgf /cnr 
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Figure  19.  Statistical  Information  for 
Compression  Loading  in  Figure  16. 
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Figure  21.  Statistical  Information  for  the  Biaxial  Loading  in 
Figure  20. 
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Figure  22.  Normalized  Constrained  Moduli,  dS^/d!^,  Versus 
Stress  Ratio,  K«o  /o,-  °7o/a\\  A11  D^tinec  Element, 

Finite  Element  and  Experimental  Results.  (a)  In  the 


Direction  of  o__,  and  (b)in  the  Direction  of 
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